% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Replication "Deconstructing the Yield Curve"
% Crump and Gospodinov (2024)
% Date: 26-JUL-2024
% Code to produce Table 7 and Table 8
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all; 
clc; 
close all;
rng(12345);
addpath('functions');
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Options
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
o.exactDataBR           = true; 
o.dataSet               = 'gsw';
o.freq                  = 'q';
o.startDateSample       = 197104; 
o.endDateSample         = 200704;
%o.endDateSample         = 201801;
o.rstarBasedOnCPI       = false;
p.N                     = 60;
p.yldMats               = [1,2,4:4:60];
p.H                     = 1;
p.CG_B_BC               = 1000;
p.BH_B_BC               = 1000;
p.B                     = 5000;
p.L                     = 3;
p.M                     = 'ROT';
p.K                     = numel(p.yldMats);
% Keep this fixed to ensure we have yoy core PCE starting in 1961Q4 for EWMA
p.startDate             = 196004; 
idxPlotCG               = 51;
idxPlotBH               = 8;
idxSaveData             = true;
idxSavePlot             = true;
nwLagh4                 = 6;

if o.exactDataBR && (o.startDateSample ~= 197104 ||  o.endDateSample > 201801)
    error('Option to use exact BR data only for sample: 197104--201801');
end

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Load data
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp(' '); disp('Loading Data...'); disp(' ');
% Load data sets
dataAll = loadBondData(p.N);
dataTag = ['dataAll.', o.dataSet, '.', o.freq];
for i = {'prc', 'yld', 'ret', 'xret', 'fwd', 'dret'}, eval([i{1}, 'Raw = ', dataTag, '.', i{1}, ';']);  end
data.mats = (1:p.N);
datesBonds = eval([dataTag, '.dates']);
data.period = eval([dataTag, '.period']);

datesAll = 100*kron((1901:2026)',ones(4,1))+kron(ones(126,1),(1:4)'); %!%
tmpIdxStart = find(datesAll==p.startDate);
tmpIdxEnd = find(datesAll==o.endDateSample);
datesAll = datesAll(tmpIdxStart:tmpIdxEnd);
TT = numel(datesAll);

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Load CPI data
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
rawCPI = readtable('../input/CPILFESL.xls', 'DataRange', 'A11');
% Lag CPI by one month to match CPo approach
rawCPI.CPILFESL = [NaN; rawCPI.CPILFESL(1:end-1)];
rawCPI.coreCPIyoy = [NaN(12,1); 100*(log(rawCPI.CPILFESL(13:end))-log(rawCPI.CPILFESL(1:end-12)))];
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Load PTR data
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
rawPTR = readtable('../input/115622-V1/data/pistar_PTR.csv');
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Load PCE data
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
rawPCE = readtable('../input/115622-V1/data/PCEPILFE.csv');
rawPCE.corePCEyoy = [NaN(12,1); 100*(log(rawPCE.PCEPILFE(13:end))-log(rawPCE.PCEPILFE(1:end-12)))];
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Load T-Bill data
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if o.exactDataBR
    rawTB3mavg = readtable('../input/115622-V1/data/TB3MS.csv');
end
rawTB3 = readtable('../input/DTB3.xls', 'Range', 'A11');
rawTB6 = readtable('../input/DTB6.xls', 'Range', 'A11');

% PTR
rawDatesPTR = 100*rawPTR.Year + rawPTR.Quarter;
data.PTR = matchDates(datesAll,rawDatesPTR,rawPTR.pistar_PTR);
% PCE
tmpIdx = find(month(rawPCE.DATE)==3, 1, 'first');
tmpData = rawPCE.corePCEyoy(tmpIdx:3:end);
rawDatesPCE = 100*year(rawPCE.DATE(tmpIdx:3:end)) + quarter(rawPCE.DATE(tmpIdx:3:end));
data.corePCEyoy = matchDates(datesAll,rawDatesPCE,tmpData);
% CPI
tmpIdx = find(month(rawCPI.observation_date)==3, 1, 'first');
tmpData = rawCPI.coreCPIyoy(tmpIdx:3:end);
rawDatesCPI = 100*year(rawCPI.observation_date(tmpIdx:3:end)) + quarter(rawCPI.observation_date(tmpIdx:3:end));
data.coreCPIyoy = matchDates(datesAll,rawDatesCPI,tmpData);
% TB3
tmpIdx = find(month(rawTB3.observation_date)==3, 1, 'first');
tmpData = rawTB3.DTB3(tmpIdx:3:end);
rawDatesTB3 = 100*year(rawTB3.observation_date(tmpIdx:3:end)) + quarter(rawTB3.observation_date(tmpIdx:3:end));
data.TB3 = matchDates(datesAll,rawDatesTB3,tmpData);
if o.exactDataBR
    tmpIdx = find(month(rawTB3mavg.DATE)==3, 1, 'first');
    tmpData = rawTB3mavg.TB3MS(tmpIdx:3:end);
    rawDatesTB3mavg = 100*year(rawTB3mavg.DATE(tmpIdx:3:end)) + quarter(rawTB3mavg.DATE(tmpIdx:3:end));
    data.TB3mavg = matchDates(datesAll,rawDatesTB3mavg,tmpData);
end
% TB6
tmpIdx = find(month(rawTB6.observation_date)==3, 1, 'first');
tmpData = rawTB6.DTB6(tmpIdx:3:end);
rawDatesTB6 = 100*year(rawTB6.observation_date(tmpIdx:3:end)) + quarter(rawTB6.observation_date(tmpIdx:3:end));
data.TB6 = matchDates(datesAll,rawDatesTB6,tmpData);
% Bond Data (replace 3-month and 6-month with Fred data)
data.yld = matchDates(datesAll,datesBonds,yldRaw(:,p.yldMats));
data.yld(:,1) = data.TB3;
data.yld(:,2) = data.TB6;

if o.exactDataBR
    warning('OFF', 'MATLAB:table:ModifiedAndSavedVarnames');    
    rawY = readtable('../input/115622-V1/data/yields.csv'); 
    yearY = floor(rawY.Date/10000);
    monthY = floor((rawY.Date-10000*floor(rawY.Date/10000))/100);
    quarterY = monthY/3;
    datesY = 100*yearY+quarterY;
    tmpY = table2array(rawY(:,2:end));
    if datesY(end)>datesAll(end)
        tmpIdx = find(datesY==datesAll(end));
        datesY = datesY(1:tmpIdx);
        tmpY = tmpY(1:tmpIdx,:);
    end
    [~,tmpIdx] = intersect(datesAll,datesY);        
    tmpYupdate = data.yld(tmpIdx,:);
    tmpYupdate(~isnan(tmpY)) = tmpY(~isnan(tmpY));
    data.yld(tmpIdx,:) = tmpYupdate;
end

[xrn,xrn4] = makeReturnsBR(data.yld);

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Make PCs
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
idxSampleStart = find(datesAll==o.startDateSample);
idxSampleEnd = find(datesAll==o.endDateSample);
T = idxSampleEnd - idxSampleStart + 1;
% Make PCs
[tmp_evecs,tmp_eigs] = eig(cov(data.yld(idxSampleStart:idxSampleEnd,:)));
[~,idx] = sort(diag(tmp_eigs), 'descend');
tmp_evecs = tmp_evecs(:,idx);
PCsLoad = tmp_evecs(:,1:3)*diag(sign(tmp_evecs(end,1:3)));
PCs = data.yld(idxSampleStart:idxSampleEnd,:)*PCsLoad;
PC1 = PCs(:,1);
PC2 = PCs(:,2);
PC3 = PCs(:,3);
% Make scaled PCs
sc = [sum(tmp_evecs(:,1)), tmp_evecs(1,2)-tmp_evecs(end,2), tmp_evecs(end,3)-2*tmp_evecs(4,3)+tmp_evecs(1,3)];
PCssc = PCs*diag(1./sc);
PC1sc = PCssc(:,1);
PC2sc = PCssc(:,2);
PC3sc = PCssc(:,3);
ehat = data.yld(idxSampleStart:idxSampleEnd,:) - PCs*PCsLoad';
s2hat = trace(ehat'*ehat)/numel(ehat);

xrnSample = xrn(idxSampleStart:idxSampleEnd,4:end);
xrn4Sample = xrn4(idxSampleStart:idxSampleEnd,4:end);
PTRSample = data.PTR(idxSampleStart:idxSampleEnd);

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Make all other variables
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if ~o.rstarBasedOnCPI
    if  o.exactDataBR
        RRQ = data.TB3mavg - data.corePCEyoy;
    else
        RRQ = data.TB3 - data.corePCEyoy;
    end
else
    RRQ = data.TB3 - data.coreCPIyoy;
end
ewmaRRQ = [NaN(4,1); ewma(RRQ(5:end),0.98)];
ewmaRRQsample = ewmaRRQ(idxSampleStart:idxSampleEnd);

CPoQ = ewma_trunc(data.coreCPIyoy,0.9615,40);
CPoQsample = CPoQ(idxSampleStart:idxSampleEnd);

datesPlot = datetime(floor(datesAll/100),4*(datesAll-100*floor(datesAll/100)),15);

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Make Figure 3(a)
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figName = figure(1);
pl0 = plot(datesPlot(idxSampleStart:idxSampleEnd),data.yld(idxSampleStart:idxSampleEnd,1), 'LineWidth', 1.25, 'Color', .7*[1,1,1]);
hold on
pl1 = plot(datesPlot(idxSampleStart:idxSampleEnd),data.yld(idxSampleStart:idxSampleEnd,2:end), 'LineWidth', 1.25, 'Color', .7*[1,1,1]);
pl3 = plot(datesPlot(idxSampleStart:idxSampleEnd),PTRSample+ewmaRRQsample, 'k', 'LineWidth', 1.25);
legend([pl0,pl3],'Yields', '$i^*$', 'Location', 'NorthEast', 'Orientation', 'Horizontal', 'Interpreter', 'Latex');
set(gca,'FontSize', 16);
hold off
if idxSavePlot && (o.endDateSample == 201801), print(figName, '-dpdf', '../output/ExampleData'); end
close all;

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Full-sample regression results
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
reg0 = nwest(mean(xrnSample(2:end,:),2),[ones(T-1,1), PCssc(1:end-1,:)],0);
reg40 = nwest(mean(xrn4Sample(5:end,:),2),[ones(T-4,1), PCssc(1:end-4,:)],nwLagh4);
reg1 = nwest(mean(xrnSample(2:end,:),2),[ones(T-1,1), PCssc(1:end-1,:) PTRSample(1:end-1)],0);
reg1_CPo = nwest(mean(xrnSample(2:end,:),2),[ones(T-1,1), PCssc(1:end-1,:) CPoQsample(1:end-1)],0);
reg41 = nwest(mean(xrn4Sample(5:end,:),2),[ones(T-4,1), PCssc(1:end-4,:) PTRSample(1:end-4)],nwLagh4); 
reg41_CPo = nwest(mean(xrn4Sample(5:end,:),2),[ones(T-4,1), PCssc(1:end-4,:) CPoQsample(1:end-4)],nwLagh4); 
reg2 = nwest(mean(xrnSample(2:end,:),2),[ones(T-1,1), PCssc(1:end-1,:) PTRSample(1:end-1) ewmaRRQsample(1:end-1)],0);
reg2_CPo = nwest(mean(xrnSample(2:end,:),2),[ones(T-1,1), PCssc(1:end-1,:) CPoQsample(1:end-1) ewmaRRQsample(1:end-1)],0);
reg42 = nwest(mean(xrn4Sample(5:end,:),2),[ones(T-4,1), PCssc(1:end-4,:) PTRSample(1:end-4) ewmaRRQsample(1:end-4)],nwLagh4); 
reg42_CPo = nwest(mean(xrn4Sample(5:end,:),2),[ones(T-4,1), PCssc(1:end-4,:) CPoQsample(1:end-4) ewmaRRQsample(1:end-4)],nwLagh4); 
reg3 = nwest(mean(xrnSample(2:end,:),2),[ones(T-1,1), PCssc(1:end-1,:) ewmaRRQsample(1:end-1)],0);
reg43 = nwest(mean(xrn4Sample(5:end,:),2),[ones(T-4,1), PCssc(1:end-4,:) ewmaRRQsample(1:end-4)],nwLagh4); 

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 1985 sample regression results
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
offset = find(datesAll==198501) - idxSampleStart;
reg0_sub = nwest(mean(xrnSample(2+offset:end,:),2),[ones(T-1-offset,1), PCssc(1+offset:end-1,:)],0);
reg40_sub = nwest(mean(xrn4Sample(5+offset:end,:),2),[ones(T-4-offset,1), PCssc(1+offset:end-4,:)],nwLagh4);
reg1_sub = nwest(mean(xrnSample(2+offset:end,:),2),[ones(T-1-offset,1), PCssc(1+offset:end-1,:) PTRSample(1+offset:end-1)],0);
reg1_sub_CPo = nwest(mean(xrnSample(2+offset:end,:),2),[ones(T-1-offset,1), PCssc(1+offset:end-1,:) CPoQsample(1+offset:end-1)],0);
reg41_sub = nwest(mean(xrn4Sample(5+offset:end,:),2),[ones(T-4-offset,1), PCssc(1+offset:end-4,:) PTRSample(1+offset:end-4)],nwLagh4); 
reg41_sub_CPo = nwest(mean(xrn4Sample(5+offset:end,:),2),[ones(T-4-offset,1), PCssc(1+offset:end-4,:) CPoQsample(1+offset:end-4)],nwLagh4); 
reg2_sub = nwest(mean(xrnSample(2+offset:end,:),2),[ones(T-1-offset,1), PCssc(1+offset:end-1,:) PTRSample(1+offset:end-1) ewmaRRQsample(1+offset:end-1)],0);
reg2_sub_CPo = nwest(mean(xrnSample(2+offset:end,:),2),[ones(T-1-offset,1), PCssc(1+offset:end-1,:) CPoQsample(1+offset:end-1) ewmaRRQsample(1+offset:end-1)],0);
reg42_sub = nwest(mean(xrn4Sample(5+offset:end,:),2),[ones(T-4-offset,1), PCssc(1+offset:end-4,:) PTRSample(1+offset:end-4) ewmaRRQsample(1+offset:end-4)],nwLagh4); 
reg42_sub_CPo = nwest(mean(xrn4Sample(5+offset:end,:),2),[ones(T-4-offset,1), PCssc(1+offset:end-4,:) CPoQsample(1+offset:end-4) ewmaRRQsample(1+offset:end-4)],nwLagh4); 
reg3_sub = nwest(mean(xrnSample(2+offset:end,:),2),[ones(T-1-offset,1), PCssc(1+offset:end-1,:) ewmaRRQsample(1+offset:end-1)],0);
reg43_sub = nwest(mean(xrn4Sample(5+offset:end,:),2),[ones(T-4-offset,1), PCssc(1+offset:end-4,:) ewmaRRQsample(1+offset:end-4)],nwLagh4); 

% p.H==1
betaHat0 = reg0.beta;
betaHat1 = reg1.beta;
betaHat1_CPo = reg1_CPo.beta;
betaHat2 = reg2.beta;
betaHat2_CPo = reg2_CPo.beta;
betaHat3 = reg3.beta;
R2_0 = reg0.rbar;
R2_1 = reg1.rbar;
R2_1_CPo = reg1_CPo.rbar;
R2_2 = reg2.rbar;
R2_2_CPo = reg2_CPo.rbar;
R2_3 = reg3.rbar;
seHat0 = reg0.beta./reg0.tstat;
seHat1 = reg1.beta./reg1.tstat;
seHat1_CPo = reg1_CPo.beta./reg1_CPo.tstat;
seHat2 = reg2.beta./reg2.tstat;
seHat2_CPo = reg2_CPo.beta./reg2_CPo.tstat;
seHat3 = reg3.beta./reg3.tstat;
tStat1 = reg1.tstat;
tStat1_CPo = reg1_CPo.tstat;
tStat2 = reg2.tstat;
tStat2_CPo = reg2_CPo.tstat;
tStat3 = reg3.tstat;
% p.H==4
betaHat40 = reg40.beta;
betaHat41 = reg41.beta;
betaHat41_CPo = reg41_CPo.beta;
betaHat42 = reg42.beta;
betaHat42_CPo = reg42_CPo.beta;
betaHat43 = reg43.beta;
R2_40 = reg40.rbar;
R2_41 = reg41.rbar;
R2_41_CPo = reg41_CPo.rbar;
R2_42 = reg42.rbar;
R2_42_CPo = reg42_CPo.rbar;
R2_43 = reg43.rbar;    
seHat40 = reg40.beta./reg40.tstat;
seHat41 = reg41.beta./reg41.tstat;
seHat41_CPo = reg41_CPo.beta./reg41_CPo.tstat;
seHat42 = reg42.beta./reg42.tstat;
seHat42_CPo = reg42_CPo.beta./reg42_CPo.tstat;
seHat43 = reg43.beta./reg43.tstat;
tStat41 = reg41.tstat;
tStat41_CPo = reg41_CPo.tstat;
tStat42 = reg42.tstat;
tStat42_CPo = reg42_CPo.tstat;
tStat43 = reg43.tstat;
% p.H==1_sub
betaHat0_sub = reg0_sub.beta;
betaHat1_sub = reg1_sub.beta;
betaHat1_sub_CPo = reg1_sub_CPo.beta;
betaHat2_sub = reg2_sub.beta;
betaHat2_sub_CPo = reg2_sub_CPo.beta;
betaHat3_sub = reg3_sub.beta;
R2_0_sub = reg0_sub.rbar;
R2_1_sub = reg1_sub.rbar;
R2_1_sub_CPo = reg1_sub_CPo.rbar;
R2_2_sub = reg2_sub.rbar;
R2_2_sub_CPo = reg2_sub_CPo.rbar;
R2_3_sub = reg3_sub.rbar;
seHat0_sub = reg0_sub.beta./reg0_sub.tstat;
seHat1_sub = reg1_sub.beta./reg1_sub.tstat;
seHat1_sub_CPo = reg1_sub_CPo.beta./reg1_sub_CPo.tstat;
seHat2_sub = reg2_sub.beta./reg2_sub.tstat;
seHat2_sub_CPo = reg2_sub_CPo.beta./reg2_sub_CPo.tstat;
seHat3_sub = reg3_sub.beta./reg3_sub.tstat;
tStat1_sub = reg1_sub.tstat;
tStat1_sub_CPo = reg1_sub_CPo.tstat;
tStat2_sub = reg2_sub.tstat;
tStat2_sub_CPo = reg2_sub_CPo.tstat;
tStat3_sub = reg3_sub.tstat;
% p.H==4_sub
betaHat40_sub = reg40_sub.beta;
betaHat41_sub = reg41_sub.beta;
betaHat41_sub_CPo = reg41_sub_CPo.beta;
betaHat42_sub = reg42_sub.beta;
betaHat42_sub_CPo = reg42_sub_CPo.beta;
betaHat43_sub = reg43_sub.beta;
R2_40_sub = reg40_sub.rbar;
R2_41_sub = reg41_sub.rbar;
R2_41_sub_CPo = reg41_sub_CPo.rbar;
R2_42_sub = reg42_sub.rbar;
R2_42_sub_CPo = reg42_sub_CPo.rbar;
R2_43_sub = reg43_sub.rbar;    
seHat40_sub = reg40_sub.beta./reg40_sub.tstat;
seHat41_sub = reg41_sub.beta./reg41_sub.tstat;
seHat41_sub_CPo = reg41_sub_CPo.beta./reg41_sub_CPo.tstat;
seHat42_sub = reg42_sub.beta./reg42_sub.tstat;
seHat42_sub_CPo = reg42_sub_CPo.beta./reg42_sub_CPo.tstat;
seHat43_sub = reg43_sub.beta./reg43_sub.tstat;
tStat41_sub = reg41_sub.tstat;
tStat41_sub_CPo = reg41_sub_CPo.tstat;
tStat42_sub = reg42_sub.tstat;
tStat42_sub_CPo = reg42_sub_CPo.tstat;
tStat43_sub = reg43_sub.tstat;

if strcmp(p.M,'ROT'), p.M=ceil((T*p.N)^(2/5)); end

if o.endDateSample == 200704
    if idxSaveData
        tmpFileName = ['BR_B',num2str(p.B) ,'_BRdata', num2str(o.exactDataBR), '_rstarCPI', num2str(o.rstarBasedOnCPI), '_', num2str(o.startDateSample),'_', num2str(o.endDateSample)];
        save(['../output/', tmpFileName, '.mat']);
    end
    return
end


% Initialize Matrices
bsBetaHat1_CG = NaN(p.B,numel(betaHat1));
bsBetaHat1_CPo_CG = NaN(p.B,numel(betaHat1));
bsBetaHat2_CG = NaN(p.B,numel(betaHat2));
bsBetaHat2_CPo_CG = NaN(p.B,numel(betaHat2));
bsBetaHat3_CG = NaN(p.B,numel(betaHat3));
bsBetaSEHat1_CG = NaN(p.B,numel(betaHat1));
bsBetaSEHat1_CPo_CG = NaN(p.B,numel(betaHat1));
bsBetaSEHat2_CG = NaN(p.B,numel(betaHat2));
bsBetaSEHat2_CPo_CG = NaN(p.B,numel(betaHat2));
bsBetaSEHat3_CG = NaN(p.B,numel(betaHat3));
bsBetaHat41_CG = NaN(p.B,numel(betaHat41));
bsBetaHat41_CPo_CG = NaN(p.B,numel(betaHat41));
bsBetaHat42_CG = NaN(p.B,numel(betaHat42));
bsBetaHat42_CPo_CG = NaN(p.B,numel(betaHat42));
bsBetaHat43_CG = NaN(p.B,numel(betaHat43));
bsBetaSEHat41_CG = NaN(p.B,numel(betaHat41));
bsBetaSEHat41_CPo_CG = NaN(p.B,numel(betaHat41));
bsBetaSEHat42_CG = NaN(p.B,numel(betaHat42));
bsBetaSEHat42_CPo_CG = NaN(p.B,numel(betaHat42));
bsBetaSEHat43_CG = NaN(p.B,numel(betaHat43));
%
bsBetaHat1_sub_CG = NaN(p.B,numel(betaHat1));
bsBetaHat1_sub_CPo_CG = NaN(p.B,numel(betaHat1));
bsBetaHat2_sub_CG = NaN(p.B,numel(betaHat2));
bsBetaHat2_sub_CPo_CG = NaN(p.B,numel(betaHat2));
bsBetaHat3_sub_CG = NaN(p.B,numel(betaHat3));
bsBetaSEHat1_sub_CG = NaN(p.B,numel(betaHat1));
bsBetaSEHat1_sub_CPo_CG = NaN(p.B,numel(betaHat1));
bsBetaSEHat2_sub_CG = NaN(p.B,numel(betaHat2));
bsBetaSEHat2_sub_CPo_CG = NaN(p.B,numel(betaHat2));
bsBetaSEHat3_sub_CG = NaN(p.B,numel(betaHat3));
bsBetaHat41_sub_CG = NaN(p.B,numel(betaHat41));
bsBetaHat41_sub_CPo_CG = NaN(p.B,numel(betaHat41));
bsBetaHat42_sub_CG = NaN(p.B,numel(betaHat42));
bsBetaHat42_sub_CPo_CG = NaN(p.B,numel(betaHat42));
bsBetaHat43_sub_CG = NaN(p.B,numel(betaHat43));
bsBetaSEHat41_sub_CG = NaN(p.B,numel(betaHat41));
bsBetaSEHat41_sub_CPo_CG = NaN(p.B,numel(betaHat41));
bsBetaSEHat42_sub_CG = NaN(p.B,numel(betaHat42));
bsBetaSEHat42_sub_CPo_CG = NaN(p.B,numel(betaHat42));
bsBetaSEHat43_sub_CG = NaN(p.B,numel(betaHat43));
%
bsR2_1minus0_CG = NaN(p.B,1);
bsR2_1_CG = NaN(p.B,1);
bsR2_1_CPo_CG = NaN(p.B,1);
bsR2_2minus1_CG = NaN(p.B,1);
bsR2_2_CG = NaN(p.B,1);
bsR2_2_CPo_CG = NaN(p.B,1);
bsR2_3_CG = NaN(p.B,1);
bsR2_41minus40_CG = NaN(p.B,1);
bsR2_41_CG = NaN(p.B,1);
bsR2_41_CPo_CG = NaN(p.B,1);
bsR2_42minus41_CG = NaN(p.B,1);
bsR2_42_CG = NaN(p.B,1);
bsR2_42_CPo_CG = NaN(p.B,1);
bsR2_43_CG = NaN(p.B,1);
%
bsBetaHat1_BH = NaN(p.B,numel(betaHat1));
bsBetaHat2_BH = NaN(p.B,numel(betaHat2));
bsBetaHat3_BH = NaN(p.B,numel(betaHat3));
bsBetaSEHat1_BH = NaN(p.B,numel(betaHat1));
bsBetaSEHat2_BH = NaN(p.B,numel(betaHat2));
bsBetaSEHat3_BH = NaN(p.B,numel(betaHat3));
bsBetaHat41_BH = NaN(p.B,numel(betaHat41));
bsBetaHat42_BH = NaN(p.B,numel(betaHat42));
bsBetaHat43_BH = NaN(p.B,numel(betaHat43));
bsBetaSEHat41_BH = NaN(p.B,numel(betaHat41));
bsBetaSEHat42_BH = NaN(p.B,numel(betaHat42));
bsBetaSEHat43_BH = NaN(p.B,numel(betaHat43));
%
bsBetaHat1_BH_CPo = NaN(p.B,numel(betaHat1));
bsBetaHat2_BH_CPo = NaN(p.B,numel(betaHat2));
bsBetaSEHat1_BH_CPo = NaN(p.B,numel(betaHat1));
bsBetaSEHat2_BH_CPo = NaN(p.B,numel(betaHat2));
bsBetaHat41_BH_CPo = NaN(p.B,numel(betaHat41));
bsBetaHat42_BH_CPo = NaN(p.B,numel(betaHat42));
bsBetaSEHat41_BH_CPo = NaN(p.B,numel(betaHat41));
bsBetaSEHat42_BH_CPo = NaN(p.B,numel(betaHat42));
%
bsBetaHat1_sub_BH = NaN(p.B,numel(betaHat1));
bsBetaHat2_sub_BH = NaN(p.B,numel(betaHat2));
bsBetaHat3_sub_BH = NaN(p.B,numel(betaHat3));
bsBetaSEHat1_sub_BH = NaN(p.B,numel(betaHat1));
bsBetaSEHat2_sub_BH = NaN(p.B,numel(betaHat2));
bsBetaSEHat3_sub_BH = NaN(p.B,numel(betaHat3));
bsBetaHat41_sub_BH = NaN(p.B,numel(betaHat41));
bsBetaHat42_sub_BH = NaN(p.B,numel(betaHat42));
bsBetaHat43_sub_BH = NaN(p.B,numel(betaHat43));
bsBetaSEHat41_sub_BH = NaN(p.B,numel(betaHat41));
bsBetaSEHat42_sub_BH = NaN(p.B,numel(betaHat42));
bsBetaSEHat43_sub_BH = NaN(p.B,numel(betaHat43));
%
bsBetaHat1_sub_BH_CPo = NaN(p.B,numel(betaHat1));
bsBetaHat2_sub_BH_CPo = NaN(p.B,numel(betaHat2));
bsBetaSEHat1_sub_BH_CPo = NaN(p.B,numel(betaHat1));
bsBetaSEHat2_sub_BH_CPo = NaN(p.B,numel(betaHat2));
bsBetaHat41_sub_BH_CPo = NaN(p.B,numel(betaHat41));
bsBetaHat42_sub_BH_CPo = NaN(p.B,numel(betaHat42));
bsBetaSEHat41_sub_BH_CPo = NaN(p.B,numel(betaHat41));
bsBetaSEHat42_sub_BH_CPo = NaN(p.B,numel(betaHat42));
%
bsR2_1minus0_sub_CG = NaN(p.B,1);
bsR2_2minus1_sub_CG = NaN(p.B,1);
bsR2_41minus40_sub_CG = NaN(p.B,1);
bsR2_42minus41_sub_CG = NaN(p.B,1);

data.fwd = NaN(size(data.yld));
data.fwd(:,1) = data.yld(:,1);
for i = 2:size(data.yld,2)
    data.fwd(:,i) = (p.yldMats(i)*data.yld(:,i) - p.yldMats(i-1)*data.yld(:,i-1))/(p.yldMats(i)-p.yldMats(i-1));
end

yldSample = data.yld(idxSampleStart:idxSampleEnd,:);
fwdSample = data.fwd(idxSampleStart:idxSampleEnd,:);

% Initial conditions for forwards (1961Q4) 
fwdSample0 = data.fwd(5,:);
% Initial conditions for yoy PCE (1961Q4) 
pred0_corePCEyoy = data.corePCEyoy(5,:);
pred0_coreCPIyoy = data.coreCPIyoy(5,:);
% Initial condition for PTR is taken from December 1961 Livingston Survey:
% Use 6m6m forward inflation at an annual rate
pred0_PTR = 100*((130.21/129.30)^2-1);
T_all = idxSampleEnd-idxSampleStart+1+40;

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Make Bauer-Hamilton bootstrap samples
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
bs1OutBH = bootstrapVAR(PCs,PTRSample,p.B,p.BH_B_BC,p.BH_B_BC,true,true);
bs2OutBH = bootstrapVAR(PCs,[PTRSample ewmaRRQsample],p.B,p.BH_B_BC,p.BH_B_BC,true,true);
bs3OutBH = bootstrapVAR(PCs,ewmaRRQsample,p.B,p.BH_B_BC,p.BH_B_BC,true,true);
bs1OutBH_CPo = bootstrapVAR(PCs,CPoQsample,p.B,p.BH_B_BC,p.BH_B_BC,true,true);
bs2OutBH_CPo = bootstrapVAR(PCs,[CPoQsample,ewmaRRQsample],p.B,p.BH_B_BC,p.BH_B_BC,true,true);
%
bs1OutBH_sub = bootstrapVAR(PCs(offset+1:end,:),PTRSample(offset+1:end,:),p.B,p.BH_B_BC,p.BH_B_BC,true,true);
bs2OutBH_sub = bootstrapVAR(PCs(offset+1:end,:),[PTRSample(offset+1:end,:) ewmaRRQsample(offset+1:end,:)],p.B,p.BH_B_BC,p.BH_B_BC,true,true);
bs3OutBH_sub = bootstrapVAR(PCs(offset+1:end,:),ewmaRRQsample(offset+1:end,:),p.B,p.BH_B_BC,p.BH_B_BC,true,true);
bs1OutBH_sub_CPo = bootstrapVAR(PCs(offset+1:end,:),CPoQsample(offset+1:end,:),p.B,p.BH_B_BC,p.BH_B_BC,true,true);
bs2OutBH_sub_CPo = bootstrapVAR(PCs(offset+1:end,:),[CPoQsample(offset+1:end,:),ewmaRRQsample(offset+1:end,:)],p.B,p.BH_B_BC,p.BH_B_BC,true,true);

for b = 1:p.B

    % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % CG Bootstrap: BR Specification (PTR, r-star)
    % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
        % BR Specification
        bsCG1 = bootstrapCG(fwdSample,PTRSample,fwdSample0,pred0_PTR,T_all,idxSampleStart-4,idxSampleEnd-4,p);
        bsCG2 = bootstrapCG(fwdSample,[PTRSample data.corePCEyoy(idxSampleStart:idxSampleEnd)],fwdSample0,[pred0_PTR,pred0_corePCEyoy],T_all,idxSampleStart-4,idxSampleEnd-4,p);
        bsCG3 = bootstrapCG(fwdSample,data.corePCEyoy(idxSampleStart:idxSampleEnd),fwdSample0,pred0_corePCEyoy,T_all,idxSampleStart-4,idxSampleEnd-4,p);
        % CPo Specification
        bsCG1_CPo = bootstrapCG(fwdSample,data.coreCPIyoy(idxSampleStart:idxSampleEnd),fwdSample0,pred0_coreCPIyoy,T_all,idxSampleStart-4,idxSampleEnd-4,p);
        if ~o.rstarBasedOnCPI
            bsCG2_CPo = bootstrapCG(fwdSample,[data.coreCPIyoy(idxSampleStart:idxSampleEnd) data.corePCEyoy(idxSampleStart:idxSampleEnd)],fwdSample0,[pred0_coreCPIyoy,pred0_corePCEyoy],T_all,idxSampleStart-4,idxSampleEnd-4,p);
        else
            bsCG2_CPo = bootstrapCG(fwdSample,data.coreCPIyoy(idxSampleStart:idxSampleEnd),fwdSample0,pred0_coreCPIyoy,T_all,idxSampleStart-4,idxSampleEnd-4,p);
        end

        % Make r-star
        bs2Rstar  = ewma(bsCG2.bsYldCG(:,1)-bsCG2.bsPred(:,2),0.98);
        bs3Rstar  = ewma(bsCG3.bsYldCG(:,1)-bsCG3.bsPred(:,1),0.98);
        bs2Rstar_CPo = ewma(bsCG2_CPo.bsYldCG(:,1)-bsCG2_CPo.bsPred(:,end),0.98);
        % Subtract 4 because start date for datesAll is 1960Q4 whereas our
        % bootstrap starts in 1961Q4 of bootstrap world
        bs2RstarSample = bs2Rstar(idxSampleStart-4:idxSampleEnd-4);
        bs3RstarSample = bs3Rstar(idxSampleStart-4:idxSampleEnd-4);
        bs2RstarSample_CPo = bs2Rstar_CPo(idxSampleStart-4:idxSampleEnd-4);
        %
        bs1PTRSample = bsCG1.bsPred(idxSampleStart-4:idxSampleEnd-4,1);
        bs2PTRSample = bsCG2.bsPred(idxSampleStart-4:idxSampleEnd-4,1);
        %
        % Make CPo Factor
        bs1CPoQ  = ewma_trunc(bsCG1_CPo.bsPred(:,1),0.9615,40);
        bs2CPoQ  = ewma_trunc(bsCG2_CPo.bsPred(:,1),0.9615,40);
        bs1CPoQSample = bs1CPoQ(idxSampleStart-4:idxSampleEnd-4,1);        
        bs2CPoQSample = bs2CPoQ(idxSampleStart-4:idxSampleEnd-4,1);                
        %
        bs1XrnSample = bsCG1.bsXrn(idxSampleStart-4:idxSampleEnd-4,4:end);
        bs1Xrn4Sample = bsCG1.bsXrn4(idxSampleStart-4:idxSampleEnd-4,4:end);
        bs2XrnSample = bsCG2.bsXrn(idxSampleStart-4:idxSampleEnd-4,4:end);
        bs2Xrn4Sample = bsCG2.bsXrn4(idxSampleStart-4:idxSampleEnd-4,4:end);
        bs3XrnSample = bsCG3.bsXrn(idxSampleStart-4:idxSampleEnd-4,4:end);
        bs3Xrn4Sample = bsCG3.bsXrn4(idxSampleStart-4:idxSampleEnd-4,4:end);
        bs1XrnSample_CPo = bsCG1_CPo.bsXrn(idxSampleStart-4:idxSampleEnd-4,4:end);
        bs1Xrn4Sample_CPo = bsCG1_CPo.bsXrn4(idxSampleStart-4:idxSampleEnd-4,4:end);
        bs2XrnSample_CPo = bsCG2_CPo.bsXrn(idxSampleStart-4:idxSampleEnd-4,4:end);
        bs2Xrn4Sample_CPo = bsCG2_CPo.bsXrn4(idxSampleStart-4:idxSampleEnd-4,4:end);
        %
        bs1PCsscSample = bsCG1.bsPCssc(idxSampleStart-4:idxSampleEnd-4,:);
        bs2PCsscSample = bsCG2.bsPCssc(idxSampleStart-4:idxSampleEnd-4,:);
        bs3PCsscSample = bsCG3.bsPCssc(idxSampleStart-4:idxSampleEnd-4,:);
        bs1PCsscSample_CPo = bsCG1_CPo.bsPCssc(idxSampleStart-4:idxSampleEnd-4,:);
        bs2PCsscSample_CPo = bsCG2_CPo.bsPCssc(idxSampleStart-4:idxSampleEnd-4,:);

        % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Make Figure 3(c)
        % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        if (b==idxPlotCG)        
            figName = figure(2);           
            pl0 = plot(bsCG2.bsYldCG(idxSampleStart-4:idxSampleEnd-4,1), 'LineWidth', 1.25, 'Color', .7*[1,1,1]);            
            hold on
            pl1 = plot(bsCG2.bsYldCG(idxSampleStart-4:idxSampleEnd-4,2:end), 'LineWidth', 1.25, 'Color', .7*[1,1,1]);            
            pl3 = plot(bs2PTRSample+bs2RstarSample, 'k','LineWidth', 1.25);
            set(gca,'FontSize', 16);
            hold off
            if idxSavePlot, print(figName, '-dpdf', '../output/ExampleCG'); end
        end

        bsReg1_0 = nwest(mean(bs1XrnSample(2:end,:),2),[ones(T-1,1), bs1PCsscSample(1:end-1,:)],0);
        bsReg41_40 = nwest(mean(bs1Xrn4Sample(5:end,:),2),[ones(T-4,1), bs1PCsscSample(1:end-4,:)],nwLagh4); 
        bsReg1 = nwest(mean(bs1XrnSample(2:end,:),2),[ones(T-1,1), bs1PCsscSample(1:end-1,:) bs1PTRSample(1:end-1)],0);
        bsReg41 = nwest(mean(bs1Xrn4Sample(5:end,:),2),[ones(T-4,1), bs1PCsscSample(1:end-4,:), bs1PTRSample(1:end-4)],nwLagh4); 
        bsReg2 = nwest(mean(bs2XrnSample(2:end,:),2),[ones(T-1,1), bs2PCsscSample(1:end-1,:), bs2PTRSample(1:end-1) bs2RstarSample(1:end-1)],0);
        bsReg42 = nwest(mean(bs2Xrn4Sample(5:end,:),2),[ones(T-4,1), bs2PCsscSample(1:end-4,:) bs2PTRSample(1:end-4) bs2RstarSample(1:end-4)],nwLagh4); 
        bsReg2_1 = nwest(mean(bs2XrnSample(2:end,:),2),[ones(T-1,1), bs2PCsscSample(1:end-1,:), bs2PTRSample(1:end-1)],0);        
        bsReg42_41 = nwest(mean(bs2Xrn4Sample(5:end,:),2),[ones(T-4,1), bs2PCsscSample(1:end-4,:) bs2PTRSample(1:end-4)],nwLagh4); 
        bsReg3 = nwest(mean(bs3XrnSample(2:end,:),2),[ones(T-1,1), bs3PCsscSample(1:end-1,:), bs3RstarSample(1:end-1)],0);
        bsReg43 = nwest(mean(bs3Xrn4Sample(5:end,:),2),[ones(T-4,1), bs3PCsscSample(1:end-4,:) bs3RstarSample(1:end-4)],nwLagh4);
        bsReg1_CPo = nwest(mean(bs1XrnSample_CPo(2:end,:),2),[ones(T-1,1), bs1PCsscSample_CPo(1:end-1,:) bs1CPoQSample(1:end-1)],0);
        bsReg41_CPo = nwest(mean(bs1Xrn4Sample_CPo(5:end,:),2),[ones(T-4,1), bs1PCsscSample_CPo(1:end-4,:), bs1CPoQSample(1:end-4)],nwLagh4); 
        bsReg2_CPo = nwest(mean(bs2XrnSample_CPo(2:end,:),2),[ones(T-1,1), bs2PCsscSample_CPo(1:end-1,:) bs2CPoQSample(1:end-1) bs2RstarSample_CPo(1:end-1)],0);
        bsReg42_CPo = nwest(mean(bs2Xrn4Sample_CPo(5:end,:),2),[ones(T-4,1), bs2PCsscSample_CPo(1:end-4,:), bs2CPoQSample(1:end-4) bs2RstarSample_CPo(1:end-4)],nwLagh4);         
        %
        bsReg1_0_sub = nwest(mean(bs1XrnSample(2+offset:end,:),2),[ones(T-1-offset,1), bs1PCsscSample(1+offset:end-1,:)],0);
        bsReg41_40_sub = nwest(mean(bs1Xrn4Sample(5+offset:end,:),2),[ones(T-4-offset,1), bs1PCsscSample(1+offset:end-4,:)],nwLagh4); 
        bsReg1_sub = nwest(mean(bs1XrnSample(2+offset:end,:),2),[ones(T-1-offset,1), bs1PCsscSample(1+offset:end-1,:) bs1PTRSample(1+offset:end-1)],0);
        bsReg41_sub = nwest(mean(bs1Xrn4Sample(5+offset:end,:),2),[ones(T-4-offset,1), bs1PCsscSample(1+offset:end-4,:), bs1PTRSample(1+offset:end-4)],nwLagh4); 
        bsReg2_sub = nwest(mean(bs2XrnSample(2+offset:end,:),2),[ones(T-1-offset,1), bs2PCsscSample(1+offset:end-1,:), bs2PTRSample(1+offset:end-1) bs2RstarSample(1+offset:end-1)],0);
        bsReg42_sub = nwest(mean(bs2Xrn4Sample(5+offset:end,:),2),[ones(T-4-offset,1), bs2PCsscSample(1+offset:end-4,:) bs2PTRSample(1+offset:end-4) bs2RstarSample(1+offset:end-4)],nwLagh4); 
        bsReg2_1_sub = nwest(mean(bs2XrnSample(2+offset:end,:),2),[ones(T-1-offset,1), bs2PCsscSample(1+offset:end-1,:), bs2PTRSample(1+offset:end-1)],0);
        bsReg42_41_sub = nwest(mean(bs2Xrn4Sample(5+offset:end,:),2),[ones(T-4-offset,1), bs2PCsscSample(1+offset:end-4,:) bs2PTRSample(1+offset:end-4)],nwLagh4);         
        bsReg3_sub = nwest(mean(bs3XrnSample(2+offset:end,:),2),[ones(T-1-offset,1), bs3PCsscSample(1+offset:end-1,:), bs3RstarSample(1+offset:end-1)],0);
        bsReg43_sub = nwest(mean(bs3Xrn4Sample(5+offset:end,:),2),[ones(T-4-offset,1), bs3PCsscSample(1+offset:end-4,:) bs3RstarSample(1+offset:end-4)],nwLagh4);
        bsReg1_sub_CPo = nwest(mean(bs1XrnSample_CPo(2+offset:end,:),2),[ones(T-1-offset,1), bs1PCsscSample_CPo(1+offset:end-1,:) bs1CPoQSample(1+offset:end-1)],0);
        bsReg41_sub_CPo = nwest(mean(bs1Xrn4Sample_CPo(5+offset:end,:),2),[ones(T-4-offset,1), bs1PCsscSample_CPo(1+offset:end-4,:), bs1CPoQSample(1+offset:end-4)],nwLagh4); 
        bsReg2_sub_CPo = nwest(mean(bs2XrnSample_CPo(2+offset:end,:),2),[ones(T-1-offset,1), bs2PCsscSample_CPo(1+offset:end-1,:) bs2CPoQSample(1+offset:end-1) bs2RstarSample_CPo(1+offset:end-1)],0);
        bsReg42_sub_CPo = nwest(mean(bs2Xrn4Sample_CPo(5+offset:end,:),2),[ones(T-4-offset,1), bs2PCsscSample_CPo(1+offset:end-4,:), bs2CPoQSample(1+offset:end-4) bs2RstarSample_CPo(1+offset:end-4)],nwLagh4);         
                
        bsBetaHat1_CG(b,:) = bsReg1.beta;
        bsBetaHat2_CG(b,:) = bsReg2.beta;
        bsBetaHat3_CG(b,:) = bsReg3.beta;
        bsBetaSEHat1_CG(b,:) = bsReg1.beta./bsReg1.tstat;
        bsBetaSEHat2_CG(b,:) = bsReg2.beta./bsReg2.tstat;
        bsBetaSEHat3_CG(b,:) = bsReg3.beta./bsReg3.tstat;
        bsBetaHat1_CPo_CG(b,:) = bsReg1_CPo.beta;
        bsBetaHat2_CPo_CG(b,:) = bsReg2_CPo.beta;
        bsBetaSEHat1_CPo_CG(b,:) = bsReg1_CPo.beta./bsReg1_CPo.tstat;
        bsBetaSEHat2_CPo_CG(b,:) = bsReg2_CPo.beta./bsReg2_CPo.tstat;        
        %
        bsBetaHat41_CG(b,:) = bsReg41.beta;
        bsBetaHat42_CG(b,:) = bsReg42.beta;
        bsBetaHat43_CG(b,:) = bsReg43.beta;
        bsBetaSEHat41_CG(b,:) = bsReg41.beta./bsReg41.tstat;
        bsBetaSEHat42_CG(b,:) = bsReg42.beta./bsReg42.tstat;
        bsBetaSEHat43_CG(b,:) = bsReg43.beta./bsReg43.tstat;
        bsBetaHat41_CPo_CG(b,:) = bsReg41_CPo.beta;
        bsBetaHat42_CPo_CG(b,:) = bsReg42_CPo.beta;
        bsBetaSEHat41_CPo_CG(b,:) = bsReg41_CPo.beta./bsReg41_CPo.tstat;
        bsBetaSEHat42_CPo_CG(b,:) = bsReg42_CPo.beta./bsReg42_CPo.tstat;        
        %
        bsR2_1minus0_CG(b) = bsReg1.rbar - bsReg1_0.rbar;
        bsR2_1_CG(b) = bsReg1.rbar;
        bsR2_2minus1_CG(b) = bsReg2.rbar - bsReg2_1.rbar;
        bsR2_2_CG(b) = bsReg2.rbar;
        bsR2_3_CG(b) = bsReg3.rbar;
        bsR2_41minus40_CG(b) = bsReg41.rbar - bsReg41_40.rbar;
        bsR2_41_CG(b) = bsReg41.rbar;
        bsR2_42minus41_CG(b) = bsReg42.rbar - bsReg42_41.rbar;
        bsR2_42_CG(b) = bsReg42.rbar;
        bsR2_43_CG(b) = bsReg43.rbar;
        bsR2_1_CPo_CG(b) = bsReg1_CPo.rbar;
        bsR2_2_CPo_CG(b) = bsReg2_CPo.rbar;
        bsR2_41_CPo_CG(b) = bsReg41_CPo.rbar;
        bsR2_42_CPo_CG(b) = bsReg42_CPo.rbar;  
        %
        bsBetaHat1_sub_CG(b,:) = bsReg1_sub.beta;
        bsBetaHat2_sub_CG(b,:) = bsReg2_sub.beta;
        bsBetaHat3_sub_CG(b,:) = bsReg3_sub.beta;
        bsBetaSEHat1_sub_CG(b,:) = bsReg1_sub.beta./bsReg1_sub.tstat;
        bsBetaSEHat2_sub_CG(b,:) = bsReg2_sub.beta./bsReg2_sub.tstat;
        bsBetaSEHat3_sub_CG(b,:) = bsReg3_sub.beta./bsReg3_sub.tstat;
        bsBetaHat1_sub_CPo_CG(b,:) = bsReg1_sub_CPo.beta;
        bsBetaHat2_sub_CPo_CG(b,:) = bsReg2_sub_CPo.beta;
        bsBetaSEHat1_sub_CPo_CG(b,:) = bsReg1_sub_CPo.beta./bsReg1_sub_CPo.tstat;
        bsBetaSEHat2_sub_CPo_CG(b,:) = bsReg2_sub_CPo.beta./bsReg2_sub_CPo.tstat;        
        %
        bsBetaHat41_sub_CG(b,:) = bsReg41_sub.beta;
        bsBetaHat42_sub_CG(b,:) = bsReg42_sub.beta;
        bsBetaHat43_sub_CG(b,:) = bsReg43_sub.beta;
        bsBetaSEHat41_sub_CG(b,:) = bsReg41_sub.beta./bsReg41_sub.tstat;
        bsBetaSEHat42_sub_CG(b,:) = bsReg42_sub.beta./bsReg42_sub.tstat;
        bsBetaSEHat43_sub_CG(b,:) = bsReg43_sub.beta./bsReg43_sub.tstat;
        bsBetaHat41_sub_CPo_CG(b,:) = bsReg41_sub_CPo.beta;
        bsBetaHat42_sub_CPo_CG(b,:) = bsReg42_sub_CPo.beta;
        bsBetaSEHat41_sub_CPo_CG(b,:) = bsReg41_sub_CPo.beta./bsReg41_sub_CPo.tstat;
        bsBetaSEHat42_sub_CPo_CG(b,:) = bsReg42_sub_CPo.beta./bsReg42_sub_CPo.tstat;  
        %
        bsR2_1minus0_sub_CG(b) = bsReg1_sub.rbar - bsReg1_0_sub.rbar;
        bsR2_2minus1_sub_CG(b) = bsReg2_sub.rbar - bsReg2_1_sub.rbar;
        bsR2_41minus40_sub_CG(b) = bsReg41_sub.rbar - bsReg41_40_sub.rbar;        
        bsR2_42minus41_sub_CG(b) = bsReg42_sub.rbar - bsReg42_41_sub.rbar;
                          
        % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%        
        % BH Bootstrap (PTR, r-star) & (CPo, r-star)
        % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%        
        bs1YldBH = bs1OutBH.bsY1(:,:,b)*PCsLoad' + sqrt(s2hat)*randn(T,numel(p.yldMats));
        bs2YldBH = bs2OutBH.bsY1(:,:,b)*PCsLoad' + sqrt(s2hat)*randn(T,numel(p.yldMats));
        bs3YldBH = bs3OutBH.bsY1(:,:,b)*PCsLoad' + sqrt(s2hat)*randn(T,numel(p.yldMats));
        bs1YldBH_CPo = bs1OutBH_CPo.bsY1(:,:,b)*PCsLoad' + sqrt(s2hat)*randn(T,numel(p.yldMats));
        bs2YldBH_CPo = bs2OutBH_CPo.bsY1(:,:,b)*PCsLoad' + sqrt(s2hat)*randn(T,numel(p.yldMats));        
        %
        bs1YldBH_sub = bs1OutBH_sub.bsY1(:,:,b)*PCsLoad' + sqrt(s2hat)*randn(T-offset,numel(p.yldMats));
        bs2YldBH_sub = bs2OutBH_sub.bsY1(:,:,b)*PCsLoad' + sqrt(s2hat)*randn(T-offset,numel(p.yldMats));
        bs3YldBH_sub = bs3OutBH_sub.bsY1(:,:,b)*PCsLoad' + sqrt(s2hat)*randn(T-offset,numel(p.yldMats));
        bs1YldBH_sub_CPo = bs1OutBH_sub_CPo.bsY1(:,:,b)*PCsLoad' + sqrt(s2hat)*randn(T-offset,numel(p.yldMats));
        bs2YldBH_sub_CPo = bs2OutBH_sub_CPo.bsY1(:,:,b)*PCsLoad' + sqrt(s2hat)*randn(T-offset,numel(p.yldMats));        

        % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Make Figure 3(b)
        % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%        
        if (b==idxPlotBH)        
            figName = figure(3);            
            pl0 = plot(bs2YldBH(:,1), 'LineWidth', 1.25, 'Color', .7*[1,1,1]);
            hold on
            pl1 = plot(bs2YldBH(:,2:end), 'LineWidth', 1.25, 'Color', .7*[1,1,1]);
            pl3 = plot(bs2OutBH.bsY2(:,1,b)+bs2OutBH.bsY2(:,2,b), 'k', 'LineWidth', 1.25);
            set(gca,'FontSize', 16);
            hold off
            if idxSavePlot, print(figName, '-dpdf', '../output/ExampleBH'); end
        end

        [bs1XrnBH,bs1Xrn4BH] = makeReturnsBR(bs1YldBH);
        [bs2XrnBH,bs2Xrn4BH] = makeReturnsBR(bs2YldBH);
        [bs3XrnBH,bs3Xrn4BH] = makeReturnsBR(bs3YldBH);
        [bs1XrnBH_CPo,bs1Xrn4BH_CPo] = makeReturnsBR(bs1YldBH_CPo);
        [bs2XrnBH_CPo,bs2Xrn4BH_CPo] = makeReturnsBR(bs2YldBH_CPo);
        %
        [bs1XrnBH_sub,bs1Xrn4BH_sub] = makeReturnsBR(bs1YldBH_sub);
        [bs2XrnBH_sub,bs2Xrn4BH_sub] = makeReturnsBR(bs2YldBH_sub);
        [bs3XrnBH_sub,bs3Xrn4BH_sub] = makeReturnsBR(bs3YldBH_sub);
        [bs1XrnBH_sub_CPo,bs1Xrn4BH_sub_CPo] = makeReturnsBR(bs1YldBH_sub_CPo);
        [bs2XrnBH_sub_CPo,bs2Xrn4BH_sub_CPo] = makeReturnsBR(bs2YldBH_sub_CPo);

        bsReg1_BH = nwest(mean(bs1XrnBH(2:end,4:end),2),[ones(T-1,1), bs1OutBH.bsY1(1:end-1,:,b) bs1OutBH.bsY2(1:end-1,:,b)],0);
        bsReg41_BH = nwest(mean(bs1Xrn4BH(5:end,4:end),2),[ones(T-4,1), bs1OutBH.bsY1(1:end-4,:,b) bs1OutBH.bsY2(1:end-4,:,b)],nwLagh4);
        bsReg2_BH = nwest(mean(bs2XrnBH(2:end,4:end),2),[ones(T-1,1), bs2OutBH.bsY1(1:end-1,:,b) bs2OutBH.bsY2(1:end-1,:,b)],0);
        bsReg42_BH = nwest(mean(bs2Xrn4BH(5:end,4:end),2),[ones(T-4,1), bs2OutBH.bsY1(1:end-4,:,b) bs2OutBH.bsY2(1:end-4,:,b)],nwLagh4);
        bsReg3_BH = nwest(mean(bs3XrnBH(2:end,4:end),2),[ones(T-1,1), bs3OutBH.bsY1(1:end-1,:,b) bs3OutBH.bsY2(1:end-1,:,b)],0);
        bsReg43_BH = nwest(mean(bs3Xrn4BH(5:end,4:end),2),[ones(T-4,1), bs3OutBH.bsY1(1:end-4,:,b) bs3OutBH.bsY2(1:end-4,:,b)],nwLagh4);
        bsReg1_BH_CPo = nwest(mean(bs1XrnBH_CPo(2:end,4:end),2),[ones(T-1,1), bs1OutBH_CPo.bsY1(1:end-1,:,b) bs1OutBH_CPo.bsY2(1:end-1,:,b)],0);
        bsReg41_BH_CPo = nwest(mean(bs1Xrn4BH_CPo(5:end,4:end),2),[ones(T-4,1), bs1OutBH_CPo.bsY1(1:end-4,:,b) bs1OutBH_CPo.bsY2(1:end-4,:,b)],nwLagh4);
        bsReg2_BH_CPo = nwest(mean(bs2XrnBH_CPo(2:end,4:end),2),[ones(T-1,1), bs2OutBH_CPo.bsY1(1:end-1,:,b) bs2OutBH_CPo.bsY2(1:end-1,:,b)],0);
        bsReg42_BH_CPo = nwest(mean(bs2Xrn4BH_CPo(5:end,4:end),2),[ones(T-4,1), bs2OutBH_CPo.bsY1(1:end-4,:,b) bs2OutBH_CPo.bsY2(1:end-4,:,b)],nwLagh4);
        %
        bsReg1_sub_BH = nwest(mean(bs1XrnBH_sub(2:end,4:end),2),[ones(T-1-offset,1), bs1OutBH_sub.bsY1(1:end-1,:,b) bs1OutBH_sub.bsY2(1:end-1,:,b)],0);
        bsReg41_sub_BH = nwest(mean(bs1Xrn4BH_sub(5:end,4:end),2),[ones(T-4-offset,1), bs1OutBH_sub.bsY1(1:end-4,:,b) bs1OutBH_sub.bsY2(1:end-4,:,b)],nwLagh4);
        bsReg2_sub_BH = nwest(mean(bs2XrnBH_sub(2:end,4:end),2),[ones(T-1-offset,1), bs2OutBH_sub.bsY1(1:end-1,:,b) bs2OutBH_sub.bsY2(1:end-1,:,b)],0);
        bsReg42_sub_BH = nwest(mean(bs2Xrn4BH_sub(5:end,4:end),2),[ones(T-4-offset,1), bs2OutBH_sub.bsY1(1:end-4,:,b) bs2OutBH_sub.bsY2(1:end-4,:,b)],nwLagh4);
        bsReg3_sub_BH = nwest(mean(bs3XrnBH_sub(2:end,4:end),2),[ones(T-1-offset,1), bs3OutBH_sub.bsY1(1:end-1,:,b) bs3OutBH_sub.bsY2(1:end-1,:,b)],0);
        bsReg43_sub_BH = nwest(mean(bs3Xrn4BH_sub(5:end,4:end),2),[ones(T-4-offset,1), bs3OutBH_sub.bsY1(1:end-4,:,b) bs3OutBH_sub.bsY2(1:end-4,:,b)],nwLagh4);
        bsReg1_sub_BH_CPo = nwest(mean(bs1XrnBH_sub_CPo(2:end,4:end),2),[ones(T-1-offset,1), bs1OutBH_sub_CPo.bsY1(1:end-1,:,b) bs1OutBH_sub_CPo.bsY2(1:end-1,:,b)],0);
        bsReg41_sub_BH_CPo = nwest(mean(bs1Xrn4BH_sub_CPo(5:end,4:end),2),[ones(T-4-offset,1), bs1OutBH_sub_CPo.bsY1(1:end-4,:,b) bs1OutBH_sub_CPo.bsY2(1:end-4,:,b)],nwLagh4);
        bsReg2_sub_BH_CPo = nwest(mean(bs2XrnBH_sub(2:end,4:end),2),[ones(T-1-offset,1), bs2OutBH_sub_CPo.bsY1(1:end-1,:,b) bs2OutBH_sub_CPo.bsY2(1:end-1,:,b)],0);
        bsReg42_sub_BH_CPo = nwest(mean(bs2Xrn4BH_sub(5:end,4:end),2),[ones(T-4-offset,1), bs2OutBH_sub_CPo.bsY1(1:end-4,:,b) bs2OutBH_sub_CPo.bsY2(1:end-4,:,b)],nwLagh4);

        bsBetaHat1_BH(b,:) = bsReg1_BH.beta;
        bsBetaHat2_BH(b,:) = bsReg2_BH.beta;
        bsBetaHat3_BH(b,:) = bsReg3_BH.beta;        
        bsBetaHat1_BH_CPo(b,:) = bsReg1_BH_CPo.beta;
        bsBetaHat2_BH_CPo(b,:) = bsReg2_BH_CPo.beta;
        %
        bsBetaHat41_BH(b,:) = bsReg41_BH.beta;
        bsBetaHat42_BH(b,:) = bsReg42_BH.beta;
        bsBetaHat43_BH(b,:) = bsReg43_BH.beta;
        bsBetaHat41_BH_CPo(b,:) = bsReg41_BH_CPo.beta;
        bsBetaHat42_BH_CPo(b,:) = bsReg42_BH_CPo.beta;
        %
        bsBetaSEHat1_BH(b,:) = bsReg1_BH.beta./bsReg1_BH.tstat;
        bsBetaSEHat2_BH(b,:) = bsReg2_BH.beta./bsReg2_BH.tstat;
        bsBetaSEHat3_BH(b,:) = bsReg3_BH.beta./bsReg3_BH.tstat;
        bsBetaSEHat1_BH_CPo(b,:) = bsReg1_BH_CPo.beta./bsReg1_BH_CPo.tstat;
        bsBetaSEHat2_BH_CPo(b,:) = bsReg2_BH_CPo.beta./bsReg2_BH_CPo.tstat;
        %
        bsBetaSEHat41_BH(b,:) = bsReg41_BH.beta./bsReg41_BH.tstat;
        bsBetaSEHat42_BH(b,:) = bsReg42_BH.beta./bsReg42_BH.tstat;
        bsBetaSEHat43_BH(b,:) = bsReg43_BH.beta./bsReg43_BH.tstat;        
        bsBetaSEHat41_BH_CPo(b,:) = bsReg41_BH_CPo.beta./bsReg41_BH_CPo.tstat;
        bsBetaSEHat42_BH_CPo(b,:) = bsReg42_BH_CPo.beta./bsReg42_BH_CPo.tstat;

        bsBetaHat1_sub_BH(b,:) = bsReg1_sub_BH.beta;
        bsBetaHat2_sub_BH(b,:) = bsReg2_sub_BH.beta;
        bsBetaHat3_sub_BH(b,:) = bsReg3_sub_BH.beta;        
        bsBetaHat1_sub_BH_CPo(b,:) = bsReg1_sub_BH_CPo.beta;
        bsBetaHat2_sub_BH_CPo(b,:) = bsReg2_sub_BH_CPo.beta;
        %
        bsBetaHat41_sub_BH(b,:) = bsReg41_sub_BH.beta;
        bsBetaHat42_sub_BH(b,:) = bsReg42_sub_BH.beta;
        bsBetaHat43_sub_BH(b,:) = bsReg43_sub_BH.beta;
        bsBetaHat41_sub_BH_CPo(b,:) = bsReg41_sub_BH_CPo.beta;
        bsBetaHat42_sub_BH_CPo(b,:) = bsReg42_sub_BH_CPo.beta;
        %
        bsBetaSEHat1_sub_BH(b,:) = bsReg1_sub_BH.beta./bsReg1_sub_BH.tstat;
        bsBetaSEHat2_sub_BH(b,:) = bsReg2_sub_BH.beta./bsReg2_sub_BH.tstat;
        bsBetaSEHat3_sub_BH(b,:) = bsReg3_sub_BH.beta./bsReg3_sub_BH.tstat;
        bsBetaSEHat1_sub_BH_CPo(b,:) = bsReg1_sub_BH_CPo.beta./bsReg1_sub_BH_CPo.tstat;
        bsBetaSEHat2_sub_BH_CPo(b,:) = bsReg2_sub_BH_CPo.beta./bsReg2_sub_BH_CPo.tstat;
        %
        bsBetaSEHat41_sub_BH(b,:) = bsReg41_sub_BH.beta./bsReg41_sub_BH.tstat;
        bsBetaSEHat42_sub_BH(b,:) = bsReg42_sub_BH.beta./bsReg42_sub_BH.tstat;
        bsBetaSEHat43_sub_BH(b,:) = bsReg43_sub_BH.beta./bsReg43_sub_BH.tstat;        
        bsBetaSEHat41_sub_BH_CPo(b,:) = bsReg41_sub_BH_CPo.beta./bsReg41_sub_BH_CPo.tstat;
        bsBetaSEHat42_sub_BH_CPo(b,:) = bsReg42_sub_BH_CPo.beta./bsReg42_sub_BH_CPo.tstat;

end

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Make CG p-values
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
pval1 = makeBsPval(bsBetaHat1_CG,betaHat1,bsBetaSEHat1_CG,tStat1,1);
pval2 = makeBsPval(bsBetaHat2_CG,betaHat2,bsBetaSEHat2_CG,tStat2,1);
pval3 = makeBsPval(bsBetaHat3_CG,betaHat3,bsBetaSEHat3_CG,tStat3,1);
pval1_CPo = makeBsPval(bsBetaHat1_CPo_CG,betaHat1_CPo,bsBetaSEHat1_CPo_CG,tStat1_CPo,1);
pval2_CPo = makeBsPval(bsBetaHat2_CPo_CG,betaHat2_CPo,bsBetaSEHat2_CPo_CG,tStat2_CPo,1);
%
pval41 = makeBsPval(bsBetaHat41_CG,betaHat41,bsBetaSEHat41_CG,tStat41,1);
pval42 = makeBsPval(bsBetaHat42_CG,betaHat42,bsBetaSEHat42_CG,tStat42,1);
pval43 = makeBsPval(bsBetaHat43_CG,betaHat43,bsBetaSEHat43_CG,tStat43,1);
pval41_CPo = makeBsPval(bsBetaHat41_CPo_CG,betaHat41_CPo,bsBetaSEHat41_CPo_CG,tStat41_CPo,1);
pval42_CPo = makeBsPval(bsBetaHat42_CPo_CG,betaHat42_CPo,bsBetaSEHat42_CPo_CG,tStat42_CPo,1);
%
pval1_sub = makeBsPval(bsBetaHat1_sub_CG,betaHat1_sub,bsBetaSEHat1_sub_CG,tStat1_sub,1);
pval2_sub = makeBsPval(bsBetaHat2_sub_CG,betaHat2_sub,bsBetaSEHat2_sub_CG,tStat2_sub,1);
pval3_sub = makeBsPval(bsBetaHat3_sub_CG,betaHat3_sub,bsBetaSEHat3_sub_CG,tStat3_sub,1);
pval1_sub_CPo = makeBsPval(bsBetaHat1_sub_CPo_CG,betaHat1_sub_CPo,bsBetaSEHat1_sub_CPo_CG,tStat1_sub_CPo,1);
pval2_sub_CPo = makeBsPval(bsBetaHat2_sub_CPo_CG,betaHat2_sub_CPo,bsBetaSEHat2_sub_CPo_CG,tStat2_sub_CPo,1);
%
pval41_sub = makeBsPval(bsBetaHat41_sub_CG,betaHat41_sub,bsBetaSEHat41_sub_CG,tStat41_sub,1);
pval42_sub = makeBsPval(bsBetaHat42_sub_CG,betaHat42_sub,bsBetaSEHat42_sub_CG,tStat42_sub,1);
pval43_sub = makeBsPval(bsBetaHat43_sub_CG,betaHat43_sub,bsBetaSEHat43_sub_CG,tStat43_sub,1);
pval41_sub_CPo = makeBsPval(bsBetaHat41_sub_CPo_CG,betaHat41_sub_CPo,bsBetaSEHat41_sub_CPo_CG,tStat41_sub_CPo,1);
pval42_sub_CPo = makeBsPval(bsBetaHat42_sub_CPo_CG,betaHat42_sub_CPo,bsBetaSEHat42_sub_CPo_CG,tStat42_sub_CPo,1);
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Make BH p-values
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
pval1_BH = makeBsPval(bsBetaHat1_BH,zeros(size(betaHat1)),bsBetaSEHat1_BH,tStat1,5);
pval2_BH = makeBsPval(bsBetaHat2_BH,zeros(size(betaHat2)),bsBetaSEHat2_BH,tStat2,5);
pval3_BH = makeBsPval(bsBetaHat3_BH,zeros(size(betaHat3)),bsBetaSEHat3_BH,tStat3,5);
pval1_BH_CPo = makeBsPval(bsBetaHat1_BH_CPo,zeros(size(betaHat1)),bsBetaSEHat1_BH_CPo,tStat1_CPo,5);
pval2_BH_CPo = makeBsPval(bsBetaHat2_BH_CPo,zeros(size(betaHat2)),bsBetaSEHat2_BH_CPo,tStat2_CPo,5);
%
pval41_BH = makeBsPval(bsBetaHat41_BH,zeros(size(betaHat41)),bsBetaSEHat41_BH,tStat41,5);
pval42_BH = makeBsPval(bsBetaHat42_BH,zeros(size(betaHat42)),bsBetaSEHat42_BH,tStat42,5);
pval43_BH = makeBsPval(bsBetaHat43_BH,zeros(size(betaHat43)),bsBetaSEHat43_BH,tStat43,5);
pval41_BH_CPo = makeBsPval(bsBetaHat41_BH_CPo,zeros(size(betaHat41)),bsBetaSEHat41_BH_CPo,tStat41_CPo,5);
pval42_BH_CPo = makeBsPval(bsBetaHat42_BH_CPo,zeros(size(betaHat42)),bsBetaSEHat42_BH_CPo,tStat42_CPo,5);
%
pval1_sub_BH = makeBsPval(bsBetaHat1_sub_BH,zeros(size(betaHat1_sub)),bsBetaSEHat1_sub_BH,tStat1_sub,5);
pval2_sub_BH = makeBsPval(bsBetaHat2_sub_BH,zeros(size(betaHat2_sub)),bsBetaSEHat2_sub_BH,tStat2_sub,5);
pval3_sub_BH = makeBsPval(bsBetaHat3_sub_BH,zeros(size(betaHat3_sub)),bsBetaSEHat3_sub_BH,tStat3_sub,5);
pval1_sub_BH_CPo = makeBsPval(bsBetaHat1_sub_BH_CPo,zeros(size(betaHat1_sub)),bsBetaSEHat1_sub_BH_CPo,tStat1_sub_CPo,5);
pval2_sub_BH_CPo = makeBsPval(bsBetaHat2_sub_BH_CPo,zeros(size(betaHat2_sub)),bsBetaSEHat2_sub_BH_CPo,tStat2_sub_CPo,5);
%
pval41_sub_BH = makeBsPval(bsBetaHat41_sub_BH,zeros(size(betaHat41_sub)),bsBetaSEHat41_sub_BH,tStat41_sub,5);
pval42_sub_BH = makeBsPval(bsBetaHat42_sub_BH,zeros(size(betaHat42_sub)),bsBetaSEHat42_sub_BH,tStat42_sub,5);
pval43_sub_BH = makeBsPval(bsBetaHat43_sub_BH,zeros(size(betaHat43_sub)),bsBetaSEHat43_sub_BH,tStat43_sub,5);
pval41_sub_BH_CPo = makeBsPval(bsBetaHat41_sub_BH_CPo,zeros(size(betaHat41_sub)),bsBetaSEHat41_sub_BH_CPo,tStat41_sub_CPo,5);
pval42_sub_BH_CPo = makeBsPval(bsBetaHat42_sub_BH_CPo,zeros(size(betaHat42_sub)),bsBetaSEHat42_sub_BH_CPo,tStat42_sub_CPo,5);

close all;

if idxSaveData
    tmpFileName = ['BR_B',num2str(p.B) ,'_BRdata', num2str(o.exactDataBR), '_rstarCPI', num2str(o.rstarBasedOnCPI), '_', num2str(o.startDateSample),'_', num2str(o.endDateSample)];
    save(['../output/', tmpFileName, '.mat']);
end

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FUNCTIONS
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function out = ewma(y,a)
        T = numel(y);
        out = NaN(T,1);
        for t = 1:T
            if t==1
                out(t) = y(t);
            else
                out(t) = out(t-1) + (1-a)*(y(t) - out(t-1));
            end
        end
end

function x=ewma_trunc(y,alpha,fobs)
    n = length(y);
    weights = alpha.^(fobs-1:-1:0)';
    weights = weights./sum(weights); % normalize to sum up to one
    gap=zeros(n,1);
    for j=fobs:n        
        gap(j) = sum(weights.*y(j-fobs+1:j));
    end
    x = gap;
end

